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Abstract 


In this paper we present a transient, fully two-phase, non-isothermal model of carbon monoxide poisoning and oxygen bleeding in the membrane 
electrode assembly of a polymer electrolyte fuel cell. The model includes a detailed description of mass, heat and charge transport, chemisorption, 
electrochemical oxidation and heterogeneous catalysis (when oxygen is introduced). Example simulation results demonstrate the ability of the 
model to qualitatively capture the fundamental features of the poisoning process and the extent of poisoning with respect to channel temperature 
and concentration. Further examples show how the multi-step kinetics can interact with other physical phenomena such as liquid-water flooding, 
particularly in the anode. Carbon monoxide pulsing is simulated to demonstrate that the complicated reaction kinetics of oxygen bleeding can 
be captured and even predicted. It is shown that variations in the channel temperature have a convoluted effect on bleeding, and that trends in 
performance on relatively short time scales can be the precise opposite of the trends observed at steady state. We incorporate a bi-functional 
mechanism for carbon monoxide oxidation on platinum-—ruthenium catalysts, demonstrating the marked reduction in the extent of poisoning, the 
effect of variations in the platinum-—ruthenium ratio and the influence of temperature. Finally, we discuss the implications of the results, extensions 


to the model and possible avenues for experimental work. 
© 2007 Published by Elsevier B.V. 
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1. Introduction 


The proton exchange membrane fuel cell (PEMFC) is an 
emerging source of clean and efficient energy, particularly for 
application in the automotive industry. In a basic description of 
the PEMFC hydrogen (H2) is electro-oxidated at the anode pro- 
ducing protons and electrons, with simultaneous oxygen (O2) 
reduction at the cathode from the protons migrating through 
the membrane and the electrons entering through an external 
circuit. The net product is simply water. The reactants are deliv- 
ered through a carbon paper, known as the gas diffusion layer 
(GDL), and reaction takes place in catalyst layers (CL) on either 
of the membrane—the catalyst in question being platinum (Pt). 
The unit comprising GDL, CL and membrane is termed the 
membrane electrode assembly (MEA). 


* Corresponding author. Tel.: +1 778 288 4292; fax: +1 604 268 6657. 
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However, there are several factors that complicate a practical 
implementation, notably flooding of the CL (restricted ingress 
of reactants), membrane drying (decreased protonic conductiv- 
ity), inefficient use of the Pt (reaction can only occur at points 
of contact between the Pt, carbon and electrolyte) and degra- 
dation, through, for example, radical attack of the membrane, 
carbon monoxide (CO) poisoning, Pt dissolution and carbon 
corrosion. These factors are in turn influenced by the proper- 
ties of the several components of the cell and the conditions 
under which the cell is operated, for example relative humidity, 
temperature and fuel composition. In order for the protons gen- 
erated at the anode to reach the cathode freely the membrane 
must remain well hydrated. The GDL must facilitate efficient 
delivery of the reactants from the channels to the CL, but, at 
the same time, must be sufficiently hydrophobic to expel any 
build-up of liquid-water, without compromising the hydration 
level of the membrane. These issues are a selection of those 
that have formed the focus of modelling studies of PEMFC. The 
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Nomenclature 
o conductivity (S m7!) 
a specific surface area of Pt (m7!) o’ surface tension (N m7!) 
aw water activity $ potential (V) 
A specific surface area of agglomerates (m7!) Q water removal constant (m7!) 
Dads desorption coefficient (mol m~?) 
C molar concentration or specific heat capacity Subscripts 
(mol m~or Jkg~! K7!) a agglomerate or anode or adsorption 
d mean pore diameter (m) ads adsorption 
D diffusion coefficient (m? s7!) b backward (reaction) 
Eads activation energy (adsorption) (J mol~!) c cathode or capillary or CCL or condensation 
Eox activation energy (oxidation) (J mol~!) CO CO 
Eo open circuit voltage (V) d dissolved or desorption 
F Faraday’s constant (C mol!) e electrolyte 
h mass transfer coefficient (s7!) g GDL or gas 
H Henry’s constant or Heaviside function H H2 
i exchange current density (A m~?) k forward (rate constants) 
IJI Leverette function i liquid 
k thermal conductivity (W m7! K7!) M free site (Pt/Ru) 
kads adsorption coefficient (m s7!) N N2 
kox electro-chemical oxidation rate (mol m7? s7!) Ox oxidation 
L thickness (m or um) O O2 
m loading (kg m~?) p pore space 
N agglomerate density (m~?) pt platinum 
p liquid pressure (Pa) T reaction 
P gas pressure (Pa) ref reference 
q surface reaction rate (mol m~? s7!) Ru ruthenium 
r interaction parameter (J mol~!) s solid or electronic 
R universal gas constant (J mol~! K7!) T total value 
Ra agglomerate radius (m) v vapour 
Ro ORR rate (mol m~? s7!) w water 
s saturation 0 reference 
S source/sink (mol m~? s~! or W m~? or A m~?) 
t time (s) Superscripts 
T temperature (K) 7 CCL 
v velocity (m s7!) CO CO 
Veell cell voltage (V) CO—O CO oxidation 
W molar mass (kg mol~!) f film 
x distance (m or um) 8 GDL 
X molar fraction H H2 
H—O H? oxidation 

Greek letters O O2 
a transfer coefficient s surface 
B symmetry factor j equilibrium value 
y diffusion rate through films (s7!) — boundary/initial value 
ô film thickness (m or um) 
A entropy (J mol~! K7!) 
€ volume fraction 
n overpotential (V) vast majority of models are steady-state, despite the obvious 
0 contact angle or surface coverage (° or —) transient nature of cell operation (performance) in automobile 
K absolute permeability (m?) applications. Degradation phenomena are also transient, impact- 
Ka and xq adsorption and desorption coefficients (m s7!) ing performance on a long time scale, such as carbon corrosion 
À membrane water content and Pt dissolution, or on a relatively short time scale, such as 
u dynamic viscosity (kg m7!s~!) CO poisoning. Moreover, several other critical issues, notably 
v fixed charge site concentration (mol m7?) freeze—start and startup—shutdown, are strictly transient. 
p density or molar area density (kg m7? kg m~?) Given its impact, water management, particularly of the 

liquid form, has received a great deal of attention, [1-9]. 
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Simplifying assumptions are regularly adopted because of the 
complexity of the system, such as one or more of: regular- 
ized liquid-water permeability, infinitely thin CL and isothermal 
conditions [5—7,10-13]. For detailed reviews of steady-state 
modelling to the year 2004 the reader is referred to [14,15]. 
Notable exceptions to the steady-state trend are [16-21], none 
of which however include liquid-water transport or an energy 
balance, which are vital to an accurate description of flooding 
and membrane hydration. More comprehensive, fully two-phase 
models are found in [22,23]. The model in [23] is able to capture 
and predict the so-called “hysteresis” phenomenon observed in 
sweep cycling (of voltage or load), as well as its behaviour with 
respect to variations in the channel conditions. 

Degradation phenomena have likewise received little atten- 
tion in the modelling literature, particularly their incorporation 
into full MEA models. Modelling of membrane degradation is 
hampered by uncertainty in the reactions driving the structural 
changes [24]. On the other hand, due to the work of Springer et 
al. [25], Baschuk and Li [26-28] and Bhatia and Wang [29], good 
models for the kinetic mechanism of CO poisoning of the anode 
have been developed. This issue is of paramount importance 
as long as reformate gas is intended as the main source of H2, 
whether on- or off-board a vehicle, until carbon-free pathways 
are developed, [30]. The concentration of CO in reformate can 
be as high as 100 ppm, which, because of its preferential adsorp- 
tion, is extremely detrimental to performance. Of the methods 
to mitigate and reverse this process the two most widely studied 
are 


1. Introduction of O% or air into the anode stream to enhance 
(provide an additional pathway for) CO oxidation on the Pt 
surfaces [31]. 

2. Alternative catalysts, for example platinum—ruthenium 
(Pt-Ru), that promote the production of oxygen-containing 
species such as hydroxyl (OH) and therefore provide an addi- 
tional pathway for oxidation of CO from the Pt surfaces 
[32]: 


In the work of Springer et al. [25], the proposed CO poisoning 
mechanism was incorporated into a steady-state, single-phase 
model of the anode catalyst layer (ACL) and GDL. Baschuk and 
Li [27] placed the same mechanism inside a one-dimensional, 
steady-state, single-phase MEA framework, later extending their 
work to include the effects of O2 bleeding using a more sim- 
plified MEA model [28]. A three-dimensional steady-state, 
single-phase phase model can be found in [33]. Bhatia and Wang 
[29] studied the transient problem by simplifying the kinetics 
proposed by Springer et al. and neglecting the transport phe- 
nomena, restricting their attention to the ACL. A similar model, 
which includes the diffusion of gas, can be found in [34]. A dif- 
ferent approach can be found in [35], in which there is assumed 
to be a competition between linearly bonded and bridge-bonded 
adsorption, in contrast to the purely linear-bonded models in 
[25-29,33,34]. This model neglects the transport problem, is 
steady-state, employs temperature-independent rate constants 
and simplifies the CO adsorption kinetics. One of the aims of 
the latter work is to investigate the role of ruthenium (Ru), when 


alloyed with Pt, in resisting the CO poisoning. A more recent 
model of Pt-Ru, [36], takes the linear-bonded approach, con- 
centrating solely on the kinetics of CO adsorption and oxidation 
in the ACL. A detailed discussion is postponed until Section 4.2. 

The kinetics of CO poisoning, on either Pt or Pt-Ru, or 
of O2 bleeding have yet to be incorporated in a two-phase 
MEA model, even under a steady-state, isothermal assumption. 
Flooding, membrane hydration, (averaged) microscopic proper- 
ties such contact angle and pore size, and the evolution of the 
system in time will inevitably interact with and influence the 
kinetic problems—as we shall demonstrate by means of several 
examples. The aim of this paper is to incorporate these kinetic 
mechanisms in the transient, two-phase, non-isothermal model 
developed in [23] (in which its predictive capability is demon- 
strated) thereby extending CO modelling for PEMFC in several 
key areas simultaneously. The following features are included. 


1. The entire MEA with explicit models of the GDL, membrane 
and both CL. 

2. A full description of the transport of charge, mass (including 
dissolved reactants and liquid-water) and thermal energy in 
the multi-phase system. 

3. The full kinetic models of Springer et al. and Baschuk and 
Li for CO poisoning and O2 bleeding [25,28], and an aug- 
mented version of the model of Enbäck and Lindbergh for 
CO kinetics on Pt-Ru [36]. 

4. The mass-diffusion limitation in the cathode catalyst layer 
(CCL) using an agglomerate model. 

5. Microstructural characterization of each component (in an 
averaged sense) and arbitrary (operating) conditions in the 
channel. 


The resulting model has significant potential as a frame- 
work for a PEMFC code that is able to predict performance 
and degradation phenomena in realistic settings. 


2. Model assumptions and equations 


The various features and assumptions that form the basis of 
the model are now listed and briefly discussed (we refer to [23] 
for full details). 


1. One-dimensional full MEA transient. The domain includes 
the entire MEA from the interface between the gas channels 
and GDL to the membrane (see Fig. 1). Each component is 
modelled explicitly. 

2. Catalyst layers. We assume that the carbon support forms 
spherical clusters, referred to as agglomerates. Surrounding 
each agglomerate is an electrolyte layer, possibly surrounded 
by a layer of liquid-water. The pores between agglomerates 
are referred to as primary pores, distinct from the smaller 
pores between the carbon particles. 

3. Reactant and gas transport. Diffusion in the primary pores of 
the CL is assumed to conform to Fick’s law. The O2 and H2 
in the primary pores of the CL reach the Pt particles on the 
agglomerate surfaces by adsorbing at the water/electrolyte 
surface and diffusing across the surrounding layers. Devia- 
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Fig. 1. Model geometry indicating the entry points of the reactants and contam- 
inant, and the transport of protons and electrons. 


tion from Henry’s law (equilibrium) provides a driving force 
for interfacial mass transfer. As in [5,22,37] the steady-state 
Darcy’s law is employed for the convective flow of gas. 

4. Proton and electron concentration. The protons are assumed 
to be transported in the form of hydronium ions, H3O*. 
For both proton and electron transport we assume electro- 
neutrality and a pseudo steady-state (see [38]). 

5. Water. Water is considered to exist in three forms: as a dis- 
solved species, as vapour and as liquid. We assume that 
the net water produced is in liquid form. Condensation and 
evaporation are modelled using the approach in [2,4,5], and 
references therein (dictated by the deviation of the local ther- 
modynamic state from equilibrium). In a similar fashion, we 
introduce phase change between vapour and dissolved water 
and between liquid and dissolved water by considering the 
deviation from an equilibrium between the phases. 

6. Temperature. A single temperature is assumed, amounting to 
an assumption of infinite rates of heat exchange between the 
phases. 

7. Kinetics. The kinetics for CO poisoning are taken from 
[25], which have also formed the basis for the work in 
[26,27,29]. This model assumes linearly bonded CO and 
Hə on Pt, with temperature, potential and coverage depen- 
dent rate constants. Kinetics for O2 bleeding are taken from 
[28] and assume that the gas-phase oxidation of CO and 
Hə is negligible. Those on Pt—Ru are an extension of the 
model in [36], which neglects the adsorption of CO on 
Ru and the Hz kinetics entirely, and are based on the bi- 
functional mechanism proposed by Watanabe and Motoo 
[39]. 


2.1. Gas phase equations 


Let s denote liquid-water (H20!9) saturation, Co, Cn, CH, 
Cco; Cco, and Cy, respectively the concentrations of O2, N2, 
H2, CO, CO2 and vapour (H2O%??) in the pore space, and Cj.e, 
i € {O, H, CO, CO2}, the concentrations of O2, H2, CO and CO2 
in the CL electrolyte and membrane. The subscript j refers to 
anode (a) and cathode (c). Mass balances yield 


a a aC; 

5p s)Ci) pa G n vC) = — Spi (1) 
a ð aC. 

—(e(1 — s)Cx) — — (oy = vcn) =0 (2) 
ot Ox Ox © 


ð 0 aC. 
ze — s)Cy) — — Dy— = UgCy = Sce + VSad (3) 
ot ox x 


where Dp,i, Dyn and Dy, i € {O, H, CO, CO2}, are the diffusion 
coefficients of O2, H2, CO, CO2, N2 and vapour in the pore 
space (see the expressions in Table 7, in which T is temperature 
and P is gas pressure). The porosity € takes the value € = ép in 
the CL and € = €g in the GDL. 

The gas velocity is assumed to follow Darcy’s law and the 
Kozeny—Carman law is used to approximate permeability 


(4) 


where « is the absolute permeability of the GDL or CL, u 
is the dynamic viscosity of the gas, the factor (1 — s)? is the 
relative permeability, d is a mean pore diameter and K is 
the Kozeny—Carman constant. The source terms are defined 
in Table 1. The quantity v = 1800 mol m~*is the fixed-charge 
site concentration of the membrane, /pe,; are volumetric mass- 
transfer coefficients from gas to electrolyte (on the gas side) 
and H;, i € {O, H, CO, CO2}, are dimensionless Henry con- 
stants. Each Ape; is approximated based on a local Sherwood 
number of 2 (for flow past a spherical particle, [40]), Apei = 
(SA Dp iSa)/(@) = O(10°), where są is the specific surface area 
of the agglomerates. The condensation/evaporation coefficient 
hpc is defined later. 


2.2. Membrane and carbon phases 


The electrolyte volume fraction, €e, has two components: one 
from the films that coat the agglomerates (ef) and the other from 
the electrolyte contained in the agglomerate interiors (el), with 
€e = ef + e. To account for the volume change due to swelling 
(assumed to impact only the film thickness) we write e£ = dy + 
0.01264, where A is the membrane water content (mol H2O/mol 


Table 1 
Sources and sinks for the gas phase Eqs. (1)—(3) 
ACL CCL GDL Meaning 
Spi hpe,i(HiCi — Cie) hpe,i(HiCi — Cie) 0 Reactant dissolution in electrolyte 
Sce —hpe(XvP — Psat) —hpe(XyP — Psat) —hpe(XyP — Psat) Condensation/evaporation 
Sad hav(Ca — C3) hay(Ca — C3) 0 Vapour/dissolved water mass transfer 
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Table 2 
Sources and sinks for the dissolved reactants, potential and dissolved-water Eqs. (6)—(8) 

ACL CCL Meaning 

S-o —4q0,ads Rolne, T, Còe)/4 O2 consumption 
StH —4aqH,ads 0 Hp consumption 
Sr,co —44C0, ads 0 CO consumption 
St,CO> aqgc0,ox + 49C-0,0x 0 CO production 
So. aF(qH,ox + 24c0,0x) FS,0 Proton source/sink 
Sos — She — She Electron source/sink 
Sar ha(Ca — Cy) ha(Ca — Ci) Liquid/dissolved water mass transfer 
Sw —44C0,ox + 49H-0, 0x —2S,0 Water production (liquid phase) 


The last term represents liquid-water production in Eq. (11). 


SO; ) and do represents the volume fraction of film without any 
swelling. The latter quantity is related to the film thickness with- 
out swelling, 5¢,9, as shown in Table 6 (see [23] for full details). 
The Pt inside the agglomerates is assumed to be inactive and the 
combined volume fraction of the carbon, Pt and small pores, €a, 
is assumed constant. The volume fraction of primary pores is 
thus ép = 1 — €a — €e. Membrane water content satisfies 


Ca 


à = ——__———— where 
1 — 0.0126C4 


dissolved water concentration(mol m~?) 
Ca = 7 (5) 


Equations for reactant concentration in the electrolyte and 
membrane are given by 


aCi 
+) = Sri + Spi (6) 


ð ð 
3 Cie) = ax (22D 


where € = ée in the CL, e = 1 in the membrane and Dye, 
i € {O, H, CO, CO2}, are the diffusion coefficients (see Table 7) 
used with Bruggeman corrections. The source terms are defined 
in Table 2; the surface reaction rates q and the volumetric 
reaction rate Ro will be defined in Section 2.6. The quantity a 
is the volumetric specific surface area of catalyst. 

Equations for potential in the electrolyte/membrane and car- 
bon, pe and ¢s, respectively, are derived from conservation of 
charge (assuming electro-neutrality and steady-state) 


ð o 0 dds 
ee (ro £) i= (a£) 48, =0 (1) 
X 


ox Ox x 


where oe and os are the protonic and electronic conductivity 
respectively (used with Bruggeman corrections) and F is Fara- 
day’s constant. € = €e in the CL and € = | in the membrane for 
Pe, while € = 1 — ég in the GDL and € = ¢, in the CL for @s. 


The source terms are defined in Table 2 and the conductivities 
in Table 7. 

The mass balance for water dissolved in the electrolyte is 
written as follows 


3 3 ICa Sh Abe 

eC 3/2p 3/2 =-§ S 

ae- ay (< tae TAU By a val 
(8) 


in which e = 1 for the membrane and € = €e for the CL, and 
Dg is the diffusion coefficient for dissolved water, given in 
Table 7(note the Bruggeman correction). The source terms Sad 
and Sq) are defined in Tables 1 and 2 respectively, and will also 
be discussed in Section 2.5. 


2.3. Energy 


The conservation of thermal energy is expressed as follows: 


a a oT 
gL + Se (scr + e(l —s)pgCgvgT — e) 


=% (9) 
k 


where p1, Pg and ps (C1, Cg and C,) are respectively the densities 
(specific heat capacities) of liquid, gas and solid, and k (oC p) is 
the volume averaged thermal conductivity (heat capacity) 


kK=kp(1 — s)e + kise + ks(1 — €), 
pCp = EspiC1 + €(1 — s)PgCg + psCs(1 — €) (10) 


where kp, ks and kı are the thermal conductivities of the pore 
space, solid (averaged) and liquid-water respectively. In the GDL 
€ = €g and in the catalyst layers € = ép. In the membrane the 
only form of heat transport is conduction. 


Table 3 
Sources and sinks for the energy Eq. (9) 
Membrane ACL CCL GDL Meaning 
Sact 0 aFaqu,ox FnceRo(ne, T, Co) 0 Activation losses 
Srey 0 —AgTaqyox AcT Rolne, T, CO.) 0 Heat of reaction 
2 2 
ade \? 3/2_ (a9; ao dibs . 
Sohm e(%) D (2) DG o)(%) d- EPa (%) Ohmic losses 


—hg Sce 


—hg; Sce 


Heat of evaporation 
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The quantity v is the liquid-water interstitial velocity and the 
source terms Sk are defined in Table 3. In these expressions — A j 
is the entropy associated with ORR (j = c) and HOR (j = a) 
and hg is the liquid-gas enthalpy change for water. 


2.4. Liquid phase 


The mass balance of liquid-water is used with Darcy’s law 
for the convective flux 


ep, ds ə (ey ks? Op 
E E RE aes 

W; ot ox (£ n) ce + Vd + ow “i H Ox 
(11) 


where the source terms are defined in Tables 1 and 2, € = €g in 
the GDL and e€ = ép in the CL. Wj, vı, pı and p are the molar 
mass, velocity, viscosity and pressure of liquid-water (with a 
relatively permeability of s$? and absolute permeability x). Since 
p = P — pe, where pe is capillary pressure, combining the equa- 
tions above yields 


ep ðs ekp ð f 3 (dpe ds ƏP 
S 
W; ot uıWı ox ds ox ox 


)) =—Sce + vSal + Sw 
(12) 


For a discussion on the form of the permeability and capillary 
pressure the reader is referred to [4,41]. In this paper we adopt the 
widely used Leverette-function found in, for example, [9]. In the 
GDL and CL the capillary pressure is defined as (respectively) 


pexo' cos68,/=F1—s), pe=0' cos, | PKs) (13) 
Kg Ke 


where IE) = 1.417£ — 2.12é? + 1.2628? is the Leverette func- 
tion, o’ is surface tension, and 6$ and Ke (6È and Kg) are the contact 
angle and absolute permeability of the CL (GDL), respectively. 


2.5. Water balance: phase change 


Condensation/evaporation is driven by the deviation from 
equilibrium, Xy P — Psat, where Psat is the saturation pressure of 
water and the first term, in which Xy is the vapour mole fraction, 
is the partial pressure of the vapour 


ep(1 — s)Xy 
ke = 
RT 


EpSP1 
hpc = ke 7 


hpc = if XP > Peat and 


if XyP< Psa (14) 


in which ke and ke are the condensation and evaporation rate 
constants, whose values are taken from [2] (see Table 7). 

In a similar fashion, the vapour-dissolved phase-change term 
in Eqs. (3) and (8) is driven by the deviation from equilibrium 
between the vapour and dissolved water, Cg — Ci, where Ci is 
the dissolved concentration at equilibrium (from [42]) 


A* = 0.3 + 10.8aw — 1602, + 14.123, 
C¥(1 + 0.0126A*) = a* (15) 


aw = XyP/ Psat is the water vapour activity. The mass-transfer 
coefficient hay is approximated from the results in [43] 


hay = Kg — s)À 
hay = Ka — s)À 


if Ca— C <0 and 
if se =o (16) 


where kg and ke are given in Table 7. 

The equilibrium membrane water content depends on the 
environment, with either relationship (15) for contact with 
vapour or A = Ay = 16.8 (or Ca = Cg) = Ay /C1 + 0.0126A/)) 
for contact with liquid. The discontinuity between the vapour- 
saturated and liquid values is known as Schroeder’s paradox. The 
mass-transfer term Sq) in Eqs. (8) and (11), and Table 2, is decom- 
posed into terms for absorption and desorption of liquid-water to 
and from the electrolyte in the CL. When the liquid-equilibrated 
water content value Ci is reached or exceeded, it is assumed 
that desorption of water from the electrolyte takes place (as liq- 
uid), the magnitude of which is driven by Cg — Cal: Adsorption 
is assumed to take place for Ca < Cail provided s > są, where 
Sx is the immobile saturation. The term hay in Table 2 therefore 
takes the form 


ha = kaesH(Ca — Ci) + KaasH(s — sx) H(—Ca + Cap) (17) 


where H(-) is the Heaviside function and kdes and kaas are the 
coefficients of desorption and absorption, which, for simplicity, 
are assumed to be constant. Their values are chosen large enough 
that Cg does not overshoot Ci significantly; see Table 7. 


2.6. Reaction kinetics 

In this section we define the reaction rates appearing in 
Table 2. The kinetics at the cathode are approximated as follows: 
oxygen reduction reaction (ORR): 


O2 + 4H* +4e7 > 2H20 (18) 


for which we use the Butler—Volmer law (in mol m~? s7!) 


airef,O 


Rolne, T, Cò e) = 
° Oe" FCret.0 


eCo (e%a Fc /RT = e teh ne/ RT) 
se 


(19) 


with exchange current density iref,o, anodic and cathodic trans- 
fer coefficients a, and œc, reference reactant molar concentration 
Cref,o, Volumetric specific surface area of catalyst (per unit vol- 
ume of catalyst layer) a and overpotential ne. The expression 
a = ayMp/L relates a to the mass specific Pt surface area (Pt 
surface area per unit mass of Pt), apt, the Pt loading, mpt, and the 
CL thickness, L. The overpotentials in the cathode and anode, 
Nc and na respectively, are defined as 


Nc = bs — be — Eo, Na = $s — Pe (20) 


where Eo is the open circuit potential, taken to be constant. We 
use an agglomerate model for the cathode in which Coe is the 
oxygen concentration at the agglomerate surfaces. It is related 
to the bulk value, Co,e, by balancing the rate of reaction with the 
rate of diffusion of reactant through an electrolyte/water film to 
the surfaces of the agglomerates. The final form of the reaction 
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rate is then 
e%aFnc/RT = e7% Fine /RT 
Ro = 4ya'Co g 


y+ a! (e%aFnc/RT Z e=% Fnc/RT) z 


1 airef,O€e 


fo (21) 
4FCret,0 


where y is a measure of the diffusion rate through the films. It 

is approximated as follows: 

A’D 
ôi 


_ "r 
Yi + Ve 
A’ = 47 (Ra + ôe)? N 


ADe, 0 
3 p| ’ 


, where 
ôe 


yı = 


(22) 


in which ôe (ôi) is the electrolyte (water) film thickness, N is 
the number of agglomerates per unit volume, A = 4 RZN is 
the specific surface area of agglomerates (we assume that all of 
the surface area is covered), with agglomerate radius Ra, and 
D; is the diffusion coefficient of O2 through liquid-water. The 
film thicknesses, ĉe and 6), are defined in Table 6(taking into 
account electrolyte swelling). The reader is referred to [23] for 
a derivation. 

It is commonly accepted that electrochemical CO oxidation 
in acidic solutions involves surface reaction between adsorbed 
CO and an oxygen-containing species [44]. On the anode we 
use the model proposed by Springer et al. in [25], based on the 
mechanism 
Hə adsorption: Hz + 2M = 2(H-M), 

H? electro-oxidation: 2(H-M) = 2Ht + 2e7 + 2M, 
CO adsorption: CO + M = CO-M, 
CO electro-oxidation: 

CO-M + H20 = M + CO2 + 2H* + 2e7 (23) 


where M represents a free Pt site. When Op is present in the 
anode, we assume that the gas-phase oxidation of CO and H3 is 
negligible in comparison to oxidation by heterogeneous catal- 
ysis. Thus, the additional reactions are assumed to follow (see 
Baschuk and Li [28]) 


O2 adsorption: O2 + 2M = 2(O-M), 
CO-M + O-M — CO» + 2M, 
O-M + 2(H-M) — H20+2M 


(24) 


CO heterogeneous oxidation: 


Hə heterogeneous oxidation: 


which provides an additional path for CO oxidation, removing it 
from the Pt surfaces at a faster rate. Equations for the evolution 
of the site coverages 6;, i € {O, H, CO, M}, are as follows (using 
notation similar to that in [28]) 


p— = 2k 


2 H ,H 
dt adsCHe9m — Daaskaa 


2 
ads OF 


S 
qH, ads 

—2kH Oy (eM Fne/ RT = e Poh hy _ 2k 6962, 

SY Ve~-Y 


qH,ox qH—O, ox 


(25) 


déco 7 
P- = kxasCco,e6M e BrOco/RT _ bok co oti Byv6co/RT 
CO, ads 
— 2k eco sinh( Fy, /2RT) — k—G6c0 (26) 
—_——_S $< LS SS’ 
CO,ox CO—0.0x 

do O 2 O 7,0 72 H-0 2 co 

P E 2kags COM ~ Dxaskads9O = kox 009% = Kx 609co 
$$$ re 
O,ads 

(27) 
Om = 1 — fu — 60 — 8co (28) 


We are assuming: H and O% (dissociative) adsorption and 
desorption based on Langmuir isotherms (qH,ads and Go,ads); 
Butler—Volmer kinetics for electro-oxidation of H2 and CO 
(qH,ox and gco,ox); CO adsorption and desorption based on a 
Frumkin isotherm (gco, ads); and Langmuir—Hinshelwood kinet- 
ics for the heterogeneous oxidation of CO and H2 (qH—0,ox and 
4cO-0,ox). The quantity p is the molar area density of catalyst 
sites. 

The Frumkin kinetics [45] for CO adsorption and desorption 
were proposed by Springer et al. [25] (see references therein) 
based on analysis of the dependence of the critical current den- 
sity on Oco. The same approach is adopted here, including the 
dependence of the forward step on 6co, with symmetry factor 
6 and interaction parameter r. The Langmuir model assumes 
that coverage does not exceed a monolayer, that there is no 
interaction between adsorbed molecules and that the apparent 
standard free energy of adsorption is coverage-independent. The 
Frumkin model assumes that this energy varies linearly with 6co 
(each additional molecule adsorbs with less ease). The interac- 
tion parameter is the slope of this variation [46]. Note that this 
dependence was neglected by Bhatia and Wang [29], as well as 
the second-order nature of the intermediate hydrogen step pro- 
posed by Springer et al. Following Baschuk and Li [28], the rate 
constants are assumed to be of the Arrhenuis type 


-EŻ ,/RT /RT 


k = kge b = bi ge “iv (29) 
The pre-exponential factors and activation energies are listed in 
Table 4 and their values will be discussed in Section 3.1. 
When Pt is replaced with Pt—Ru, these equations require mod- 
ification. The kinetics in that case, including the adsorption of 
water and formation of hydroxyl] radicals (OH), are presented 


and discussed in Section 4.2. 


2.7. Initial-boundary conditions, numerical details and 
parameters 


For the discussion on the boundary conditions we recall 
Fig. 1. At the interfaces between the membrane and CL, x = x2 
and x = x3, the gas-phase and liquid-water fluxes are taken to 
be zero; that is, the gas species and liquid-water are assumed not 
to penetrate the membrane. Similarly, we assume that the fluxes 
of protons and dissolved species at the interfaces between the 


8 


Table 4 
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Pre-exponential factors and activation energies for the rate coefficients (29) in the system (25)—(28) 


Adsorption Desorption Electro-oxidation 
H2 kiso =3ms7! bho = 4.18 x 10!! mol m~? ko 9 = 23.1 molm~ s7! 
Ends = 10.4kJ mol™ Eris, =87.9 kJ mol”! Ek = 16.7molm~* s7! 
co KR 9 =3 x 10¢ ms"! DER 9 = 6.87 x 10° mol m~? Key = 3.4 x 10° mol m~? s7! 
EQ, = 47.3 kJ mol! ESQ, = 100 kJ mol! EO, = 127 mol m~? s7! 
O2 kess0 =6x 10" ms" Devi = 1.36 x 10° mol m~? 
EQ x = 14.2kI mol"! EO.» = 250 kJ mol! 
Heterogeneous oxidation 
H—O kH-O = 3.28 x 10f mol m~? s7! 
ERO = 65.9 kJ mol"! 
C—O kO = 8.3 x 10’ molm~ s7! 
ECO = 90kJ mol™! 
CL and GDL, x = x; and x = x4, are negligibly small given by the formula in [47] 
Iie Ibe ICa  Sdoe Ihe logig Psat = —2.1794 + 0.02953AT — 9.1837 x 10-° AT? 
io ne ae a R 41.4454 x 10-7AT3 (32) 
oC; ðs ; ; : 
nan 0 (x = x2, x3) (30) where AT = T — 273. From these relationships we obtain 


At the interfaces between the channels and GDL, the gas 
concentrations, temperature, water activity and pressure are pre- 
scribed. The cell voltage, Veen, is prescribed at the cathode 
channel/GDL interface, and at the anode channel/GDL interface 
we assume a zero concentration of electrons. 


Ci(xo) = Cre, Ci (x5) = Cia iE {H, O, CO, CO», N, v}, 
Qs(x0) = Ven, Ps(x5)=0, Txo)= Te, T&s)= T, 
P(xo) = Pe, P(x5) = Pa, aw(xo) = aw,cs dw (x5) = awa 


(31) 


The values of C v,c and Č v,a are calculated from the water activity. 
The saturation pressure in atm, a function of temperature, is 


Table 5 
Channel conditions assumed in the calculations, unless otherwise specified 


Cyc = Awe Psat,e/ RTe, where Psat,c is the cathode-channel satu- 

ration pressure. A similar calculation applies on the anode side. 
The final boundary conditions are those for liquid-water at the 

interfaces between the gas channels and the GDL. We approx- 

imate them with the following steady-state flux conditions at 

x = Xo and x = x5 (see [23] for details) 

Os 


ox 
where &2 ; = 0 corresponds to zero water removal from the anode 
channel, j = a, or cathode channel, j = c. Q j is likely to depend 
sensitively on the flow rate in the channel [48]. 

Initial conditions for pressure, temperature and vapour con- 
centrations are consistent with the conditions in the channels. 
The water content of the membrane/electrolyte is assumed to be 
given by equilibrium with vapour in the channels. The saturation 
at the initial time is assumed zero. 


-Qjs=0, jef{ac} (33) 


Symbol Quantity Size 

Te Cathode channel temperature 70°C 

Ta Anode channel temperature 70°C 

awe Cathode channel water activity 1 

awa Anode channel water activity 1 

Cove Oxygen concentration in cathode channel 18.13 mol m~? 
Cua Hydrogen concentration in anode channel 34.52 mol m~? 
P: Gas pressure in the cathode channel 300 kPa 

Pa Gas pressure in the anode channel 300 kPa 

Cyc Vapour concentration in cathode channel 15.91 mol m~? 
Cia Vapour concentration in the anode channel 15.91 mol m~? 


0.075 m~! (assumed) 


i Liquid-water removal constants (anode and cathode) 
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Table 6 
Unless otherwise specified these parameter values, relating to structural and electrochemical properties, were used in the calculations 
Symbol Quantity Size Ref. 
Structural 
L Catalyst layer thickness 25 um assumed 
Lm Membrane thickness 50 um Assumed 
Lg GDL thickness 200 um Assumed 
é Electrolyte volume fraction in agglomerates 0.15 Assumed 
Ea Volume fraction of agglomerates 0.4215 [55] 
Es Volume fraction of carbon in CL 0.2 Estimated 
Eg Porosity of the GDL 0.74 [56] 
Ra Agglomerate radius 0.5 um [55] 
5e,0 Electrolyte film thickness w/o swelling 0.1 um [23] 
dy Electrolyte volume fraction w/o swelling NTR, + 8e,0)° — R3] [23] 
ôe Electrolyte film thickness (m) i RÈ + 36 — Ra [23] 
ôi Water film thickness (m): Rs = Ra + ôe P/R + PL — Rs [23] 
N Agglomerate density 5.8 x 10!7 m~? Estimated 
apt Specific surface area of Pt 1000 cm? (mg Pt)! [57] 
Mpt Pt loading 0.4 (mg Pt) cm~? [3] 
(A (6$) Catalyst-layer (GDL) contact angle 90° (120°) [58] 
dg GDL pore size 10 um Assumed 
de Catalyst-layer pore size 2 um Assumed 
Electrochemical 
iref,O Cathode exchange current density 10-2 A m~? [22] 
Cref,0 Reference Oz concentration 0.05 mol m~? Assumed 
Qe Cathodic transfer coefficient 0.55 Measured 
Qa Anodic transfer coefficient 0.45 Measured 
B Symmetry factor 0.1 [28] 
r Interaction parameter 39.77 kJ mol`! [28] 
Ba Symmetry factor (H20) 0.5 [36] 
r2 Interaction parameter (H20) 1.34 kJ mol”! [36] 
Eo Open circuit potential 0.95 V Assumed 
p Molar area density of Pt sites 0.01042 mol m~? [29] 
PRu Molar area density of Ru sites p/2mol m~? Variable 
ky Forward rate constant (CO oxidation) 1077 mol m~? s7! Fitted 
kı Backward rate constant (CO oxidation) 107° mol m~? Fitted 
k2 Forward rate constant (H20 adsorption) 4 x 1077 mol m~? s7! [36] 
k2 Backward rate constant (H20 adsorption) 9 x 1078 mol m~? [36] 
— Ac Entropy associated with ORR 163.7 Jmol! K7! [59] 
— Aa Entropy associated with HOR? 0Jmol~!K! [59] 


The initial-boundary value problem was solved in COMSOL 
3.2a, on a uniform grid (typically 128 points) using quartic 
Lagrange polynomials as trial and test functions. The switch 
functions were substituted with hyperbolic tangent functions to 
smooth the discontinuities, a standard procedure. 

The default set of parameter values is given in Tables 5-7 . 
Parameter fitting of some of the electrochemical rate constants 
was found to be necessary, as in [25—29,33,34]. Several other 
parameters are estimated (as indicated in Tables 5-7), and the 
rest are found from the literature (with references provided). 


3. Results and discussion 
3.1. Validation and dependence on concentration 

As a confirmation that the kinetics have been correctly cap- 
tured we first compare our results to the experimental data 


collected by Bhatia and Wang [29], for temporal variations in 
the current density at a fixed cell voltage, under several CO 


concentrations. To reproduce those results we use pressures 
of 3 atm and temperatures of 80° C on both anode and cath- 
ode, with fully saturated conditions (aw = 1 in both channels). 
Operation at 0.6 V for 2h with zero CO concentration was sim- 
ulated. On completion CO was introduced in the anode stream 
at different concentrations and calculations were continued to 
simulate 2h of poisoning, for both 40%H2/N2/CO and H2/CO. 
The results are demonstrated in Fig. 2, in which it can be seen 
that for 40%H2/N2/CO the current density drops substantially, 
by around 60% for CO concentrations of 10 ppm and by around 
90% at 100 ppm. With H»/CO the degradation is less severe at 
fixed concentrations of CO. This figure is to be compared with 
the experimental results in Fig. 6 of [29]; the trends are captured 
correctly, with respect to the concentration of both CO and H2. 

Fig. 3 shows the evolution of CO surface coverage, 9co, in the 
two cases of 10 and 100 ppm in 40%H?/N2/CO. Even at 10 ppm 
the level of CO is extremely high, indicating that CO coverage 
must be close to a monolayer before its effect on the the current 
density is strongly felt, the reason for this being the (relatively) 
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Table 7 

Heat, charge and mass transfer/transport properties used in the calcu lations, unless otherwise specified 

Symbol Quantity Size Ref. 

Dp,o O72 diffusion coefficient 4.23 x 10-%(e(1 — s)T)?/? /P m°? s7! 40] 

Dow Hp diffusion coefficient 1.43 x 10-8(e(1 — s)TP?/P m? s7! 40] 

Dy Vapour diffusion coefficient 5.27 x 107° (e — s)T\/?/P m? s7! 40] 
Dp.co CO diffusion coefficient 4.12 x 10-%(e(1 — s)TY!?/P m? s~! 40] 
Dp,co, CO, diffusion coefficient 3.6 x 10-%(e(1 — TY? /P m? s7! 40] 

Dy No diffusion coefficient 3.6 x 10-%(e(1 — s)T)/? /P m? s7! 40] 

Da Diffusion coefficient for HO“ 4.17 x 1078A(1 + 161 e7?) e7?®6/T m? s7! 60] 

D; O% diffusivity in H20"4 (60°C) 4.88 x 107°? m? s7! 40] 

Deo O% diffusivity in Nafion 3.1 x 1077e72768/T m? s7! 61] 

Dey ê H diffusivity in Nafion 6.92 x 107°? m? s7! 62] 
Deco ê CO diffusivity in Nafion 4.38 x 107°? m? s7! 62] 
De,co CO, diffusivity in Nafion 4.12 x 107°? m? s7! 62] 

Npe,i Mass transfer rates 10° s7! Estimated 
Ho? O2 Henry’s law constant 0.15 16] 

Hy? Hz Henry’s law constant 0.63 6] 

Heo CO Henry’s law constant 0.32 

Hco, * CO> Henry’s law constant 1.94 

Kj (Kg) Absolute permeability of CCL (GDL) 107! (8.7 x 107!) m? 56,63] 
My H20"9 viscosity 107° Pas 4] 

H Dynamic viscosity air 1.8 x 1075 Pas Estimated 
H Dynamic viscosity H2 8.4 x 1076 Pas Estimated 
o' Surface tension 0.07 N m7! 4] 

k Catalyst layer thermal conductivity 0.67 W m7!K7! 64] 

km Membrane thermal conductivity 0.67 W m7!K7! 64] 

kg GDL thermal conductivity 1.67Wm7!K7! 64] 

pC H20" heat capacitance 4.187 x 10°Jm73K7! Estimated 
PgCg Air heat capacitance 10? Jm~°KT! Estimated 
PmCm Membrane heat capacitance 2.18 x 106 Jm~°KT! Estimated 
PeCe Carbon heat capacitance 1.61 x 106 Jm~?KT! Estimated 
—Ac Entropy associated with ORR 163.7 J mol`! K7! Lamp 
—Aa Entropy associated with HOR? OJ mol~!K~! 59] 

Os Electronic conductivity 500S m7! 65] 

Oe Protonic conductivity e!286(1/303-1/7)(9. 5142, — 0.326) S m7! 47] 

kaes H20"4 Desorption coefficient 100 Assumed 
Kads H2 0“is adsorption coefficient 10 Assumed 
Ka Absorption constant 10-6 ms! 43] 

Ka desorption constant 3.3 x 1076 ms7! 43] 

ke Evaporation coefficient 100 s~! atm7! 2] 

ke Condensation coefficient 100s~! 2] 


è Approximated by value in/for liquid-water at 60°C. 
b For a detailed discussion see [66]. 
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Fig. 2. Simulated evolution in the current density for several concentrations of CO in 40%H2/N2 and Hp. For a list of parameter values see Tables 5-7. 
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Fig. 3. Simulated profiles of CO surface coverage, 0co, for 100 ppm (left-hand figure) and 10 ppm (right-hand figure) of CO in 40%H2/N2, corresponding to the 


calculations in Fig. 2. 


rapid rate of H oxidation. The quantitative values are dependent 
on the parameters — those determining the rates of adsorption, 
desorption and oxidation — as was discussed in detail by Springer 
et al. [47]. Bhatia and Wang [29] adjusted the two parameters 
(assumed constant in their work) ke, and beO in order to match 
the experimental results. To achieve the correct qualitative trends 


a similar procedure was applied here with adjustments in kA, 0 


and kai: the latter to a value close to that in [47] where it was 
approximately four orders of magnitude smaller than in the base 
case of [28]. In fact the corresponding term (CO heterogeneous 
oxidation) was neglected entirely by Bhatia and Wang [29] since 
it was deemed too small to contribute, at least at the cell voltages 


considered. The remaining values in Table 4 are taken from [28]. 


3.2. Effects of cell voltage and sweep rate: potentiostatic 
sweep 


One of the effects not captured in the previous result (Fig. 2) 
is that of cell voltage on the extent of degradation. Of course, 
under steady-state conditions this effect is naturally contained 


cell voltage [V] 


0 0.25 0.5 


0.75 1 1.25 1.5 


current density [Acm~7] 


in the polarization curves, which represent, generally speaking, 
the behaviour at time t = 00 (i.e., sufficiently long for a ‘steady- 
state’ to be reached). In this section we present simulations of 
transient variations in the cell voltage (potentiostatic sweeps) for 
time scales smaller than those required for relaxation to steady- 
state. These experiments can yield results that are qualitatively 
similar to the steady-state curves but in general contain a richer 
structure [22,23,49]. 

The first result is shown on the left-hand side of Fig. 4, 
demonstrating the variation in performance with 40%H»/N2 and 
different CO concentrations. These curves are generated by a 
simulated voltage sweep from 1 to OV in 80min. As with the 
steady-state polarization curves (see for example [25]) there are 
two characteristic turning points for large CO concentrations, 
one occurring at high cell voltage and marking the onset of 
poorer performance compared to 0 ppm, and one at a lower cell 
voltage, at which performance improves markedly. The expla- 
nation for these turning points is apparent from Fig. 5, which 
shows the simulated profiles of CO and H2 surface coverage, 
Oco and 6y, and simulated electro-oxidation rates, agco,ox and 
aqu,ox, for 50 ppm. 
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Fig. 4. The left-hand figure shows polarization curves from simulated voltage sweeps (1-0 V in 80 min) for a range of CO concentrations in 40%H2/N2. The 
right-hand figure demonstrates the effect of the sweep rate. The remaining parameter values are given in Tables 5-7. 
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Fig. 5. The first two rows show simulated profiles of the surface coverage, 04 and @co, and the electro-oxidation rates, aqy,ox and agco,ox- The bottom row shows 
the evolutions of anode overpotential, na, and Hz surface coverage, Oy, at the interface between the ACL and membrane, x = x3. Note that these plots correspond to 
the calculation for 25 ppm CO in 40%Hp2/N>2 on the left-hand side of Fig. 4: a potential sweep from 1 to 0 V in 80 min. 


The first turning point is caused by the rapid coverage of 
sites, the extent of which increases with increasing CO concen- 
tration at fixed cell voltage (see Fig. 4). In the range of cell 
voltage between 1 V and the first turning point the rates of H2 
and CO oxidation are slow. At a cell voltage of approximately 
0.45 V the rate of H2 oxidation increases rapidly because of its 
exponential dependence on anode overpotential, the maximum 
in which evolves as shown in the bottom right of Fig. 5. This 


rapid increase leads to the second turning point in the curve in 
Fig. 4. Note that in [26] the authors postulate that increased CO 
electro-oxidation is the cause of the second turning point. This 
however seems an unlikely explanation given the magnitude of 
CO oxidation in comparison with that of H2 (see Fig. 5); larger 
values of KO will shift the location of the turning point to higher 
voltage values because of decreased CO coverage, but will not 
be the cause (see the steady-state results and discussion in [47]). 
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Fig. 6. Simulated profiles of CO surface coverage, co, and H2 electro-oxidation rate, aq4,ox, for 50 ppm in 40%H2/N2, corresponding to the calculations on the 
right-hand side of Fig. 4. The two cases shown here are for sweeps from 1 to 0 V in 10 min (left-hand column) and 80 min (right-hand column). 


Unsurprisingly, the results contained in the left-hand side of Clearly, the shorter the time the better the performance at fixed 


Fig. 4 are dependant on the rate of decrease of the cell volt- cell voltage; a manifestation of the decreasing time available to 
age (sweep rate). The right-hand side of Fig. 4 demonstrates reach CO coverage levels that inhibit the H2 oxidation substan- 
the trend as the sweep time is decreased from 80 to 10min. tially (see Fig. 6). The first turning point therefore occurs at a 
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Fig. 7. Polarization curves generated from simulated potential sweeps (1-0 V in 80 min) with 50 ppm CO in 40%H2/N2, for a range of anode-channel pressure, P,. 
Note than in the left-hand figure ay,, = dw, = 1, whereas in the right-hand figure ay, = aw, = 0.4, in which case anode flooding is not apparent. The remaining 
parameter values are given in Tables 5-7. 
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lower cell voltage as the sweep rate is decreased. This has the 
natural consequence of shortening the distance in cell-voltage 
space between the the two turning points; by the time the CO 
coverage is significant the cell voltage is low enough to bring 
about the onset of the second turning point. At a low enough 
sweep rate the extent of CO coverage is too low to engen- 
der the first turning point, as illustrated in the 10min sweep 
in Fig. 4. 
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3.3. Effects of anode channel pressure, temperature and 
humidity 


The effect of anode channel pressure on the CO poisoning 
effect was investigated experimentally by Springer in [47]. It 
was found that increasing the anode pressure had a detrimental 
effect on performance, measured at steady-state, indicating that 
the kinetic order in free catalyst sites in Eq. (25) exceeds unity 
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Fig. 8. The top row of figures shows simulated saturation profiles generated by a voltage sweep from 1 to 0 V in 80 min, with 50 ppm CO in 40%Hp/N2, for two 
different anode-channel pressures (corresponding to Fig. 7). The middle row of figures shows the profiles of Hz surface coverage and electrolyte concentration, 
Oy and Cy respectively, for the case P, = 1.5 atm, in which the anode flooding severely restricts the access of H2 to the Pt surfaces. The bottom row of figures 
demonstrates that the flooding is caused by saturation of the gas phase with vapour, with low levels of Hz concentration and aw, = 1. 
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Fig. 9. Polarization curves generated from simulated potential sweeps (1—0 V in 80 min) with 50 ppm CO in 40%H2/Nz2, for a range of anode channel humidity (left 
hand figure with T = 70° C in both channels) and a range of channel temperature (right-hand figure with ay, = aw, = 1). The remaining parameter values are 


given in Tables 5-7. 


(here assumed to be 2). For the transient case, while the same 
remains generally true for high anode pressures, the simulations 
presented next indicate that variations in anode pressure have 
a more convoluted effect, depending, in particular, on the the 
water activity in the channels. 

Fig. 7 demonstrates that the steady-state experimental results 
in [47] (Fig. 8 of that paper) are qualitatively similar to potential 
sweeps from 1 to 0 V in 80 min, particularly for cell voltage away 
from the lower limit. Deviations from the steady-state trend arise 
at low values of anode channel pressure, for low cell voltage, 
where the performance deteriorates sharply. For an explanation 
we refer to Fig. 8, the top row of which shows the evolutions of 
the saturation profiles for the two extreme cases in Fig. 7: Pa = 4 
and 1.5 atm. Clearly, the anode floods as the channel pressure is 
reduced, restricting the access of H2 to the Pt particle surfaces 
and causing the sharp drop in performance observed in Fig. 7. 
The origin of the flooding is the low concentration of H2, which 
leads to almost complete consumption in a region (in both the 
electrolyte and gas phases) close to the membrane/ACL interface 
(see the middle row of figures), locally forcing the gas pressure 
to decrease and molar fraction of vapour, Xy, to increase; the net 
result is an increasing overpressure Xy P — Psat (see the bottom 
row of Fig. 8) and therefore an increased rate of condensation, 
Sce, defined in Table | and Eq. (14). Another important effect in 
this result will be the restricted back-diffusion of water through 
the membrane (from anode to cathode), compared to the case 
of Pa = 4atm, because of the drop in the current density. This 
will lead to a greater retention of dissolved, and therefore liquid, 
water in the anode. Of course, this phenomenon is exacerbated 
by the fully humidified channel condition, awa = 1. In fact, for 
áw, = Awe = 0.4, the flooding is not seen, as demonstrated on 
the right-hand side of Fig. 7. What is also noticeable is the bet- 
ter performance at aw, = aw, = 1, which can be attributed to 
the greater membrane water content (and therefore increased 
protonic conductivity) on the cathode side. 

For high anode channel pressure the best performance is gen- 
erally seen at dw, = dw = 1, again because of the increased 
protonic conductivity, particularly on the cathode side. This is 


demonstrated in the polarization curves on the left-hand side of 
Fig. 9, for a range of channel humidity at the base-case anode 
pressure of 3 atm. Mass-transport limitations, arising from the 
cathode side, are noticeable at low cell voltage, manifested as 
an increasingly negative slope of the polarization curve: notice 
that the two cases awa = aw, = 1 and aw, = dw, = 0.8 even- 
tually cross. The results in Figs. 7 and 9 suggest that when CO is 
present in the anode stream the best performance is obtained for 
low anode pressure and high water activity in the cathode. When 
the cell is operated fully humidified in the anode channel there 
exists an optimal (low) anode pressure between 1.5 and 3 atm for 
which flooding in the anode can be avoided. Conversely, if the 
intention is to operate at an (low) anode pressure of 1.5 atm, there 
exist optimal channel water activities to avoid both flooding in 
the anode and severe mass-transport limitations in the cathode. 

A further non-trivial observation is that mass transport limi- 
tations in the cathode, typically a major concern, seem to be less 
apparent when CO is present in the anode stream, even at high 
water activity. The lower current density with a CO-fed anode is 
the likely explanation; it leads to reduced water production and, 
simultaneously, reduced back diffusion of dissolved water from 
anode to cathode. 

The effect of changes in the channel temperatures, at a fixed 
channel pressure, is shown on the right-hand side of Fig. 9. 
In general an increased temperature will enhance performance, 
primarily because of the increased membrane conductivity (see 
the form of de in Table 7) and the resultant increased overpo- 
tential; the reaction rate (19) has an exponential dependence 
on overpotential, which outstrips the reduced concentration of 
reactants in the channels, on which the reaction rates have an 
algebraic dependence. This result on its own is not surprising 
but the transient aspect of the results does however reveal an 
important detail: as the temperature is lowered, in the upper part 
of the curve better performance is attained at fixed cell voltage. 
We refer to Fig. 10, which compares the evolutions of surface 
coverage, 0y and Oco, and volumetric H? oxidation rate, aqy ox, 
at the two temperatures, evaluated at the interface between the 
membrane and ACL, x = x3. 
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Fig. 10. The evolutions of surface coverages, 64 and 6co, and H2 oxidation rate, aqH,ox, evaluated at the membrane/ACL interface, x = x3, in the two cases 
Ta = 60 and 80° C on the right-hand side of Fig. 9. These correspond to a simulated potential sweep from 1 to 0 V in 80 min, with 50 ppm CO in 40%H2/N2 and 


awa = Awe = 1. 


At the lower temperature CO (H2) coverage is lower (greater) 
up to t © 35 min, with the opposite being true beyond this time. 
This clearly influences the relative shapes of the profiles in 
Fig. 10. However, this highly non-linear behaviour is not easy 
to predict. The lower temperature leads to reduced CO and H2 
adsorption rate constants through the Arrhenius activation ener- 
gies (29) and the Frumkin term in (26). The rate constant for 
CO adsorption, ke9, experiences a more dramatic decrease than 
that for H2 adsorption, kis and therefore dominates, starting 
from the initial conditions of 6co = 04 = 0 and zero desorp- 
tion. The result is a lower CO coverage and (therefore) higher 
H coverage. The desorption rate of CO increases as 0co — 1. 
However, in this limit the ratio of desorption to adsorption is pro- 
portional to exp((—100 + r@co)/RT) ~ exp(—60/RT), which 
is much smaller at the lower temperature. In other terms, the 
long-time value of CO coverage will be greater for the lower tem- 
perature, although in the early stages the CO coverage is lower. 
In summary, there is a sensitive balance between the decreas- 
ing rate of adsorption, decreasing rate of desorption and the 
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decreased membrane hydration at the lower temperature. The 
evolution of reaction rate seen on the right-hand side of Fig. 10 
(and therefore the curves seen in Fig. 9) is determined by this 
balance in a non-trivial way. Note that the final changes in the 
slopes of these curves at around tf = 50 min coincide with the 
second turning points (H2 oxidation) in the polarization curves. 


4. Mitigation techniques 


4.1. O2 bleeding: CO pulse simulation and temperature 


effect 


A known method for preventing the degradation seen in Fig. 2 
is to inject air or neat O% into the anode fuel stream, the so- 
called O2 bleeding method. The results of a series of transient 
experiments investigating O2 bleeding can be found in [50,51], 
in which the cell voltage is measured for a constant load with 
time-varying compositions of H2, CO and O2. It was also found 
in these and other studies that a degraded cell has the capacity 
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Fig. 11. The left-hand figure demonstrates the transient effect of CO pulsing and the ability of the system to recover for the indicated levels of CO at three cell voltage 
values. The right-hand figure demonstrates the effect of O2 bleeding, at different concentrations in 40%H2/N2/CO, for a cell voltage of 0.6 V (remaining parameter 
values are given in Tables 5-7). The evolution of the surface coverages for 4% and 0.1% Oz bleeding are demonstrated in Fig. 12. 
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Fig. 12. Plots of O2, CO and H3 surface coverages evaluated at the interface between the ACL and membrane, x = x3, corresponding to the examples on the 


right-hand side of Fig. 11 for 4% and 0.1% Oz bleeding at 0.6 V. 


to recover when fed with CO-free H2, the extent of recovery 
depending on the initial input of CO and the poisoning and 
recovery times. We begin this section with simulations akin to 
the experiments in [50,51], instead calculating the current den- 
sity as the cell voltage is varied with time-varying compositions 
of H2, CO, O2 and N2. The main objective here is to demon- 
strate that qualitatively the outcome of such experiments can be 
predicted with the model. 

We simulated 1h of operation with 40%H2/N2/0 ppm CO 
in the anode and a fixed cell voltage (other parameters as 
in Tables 5-7), long enough to reach a steady current den- 
sity. A concentration of 50 ppm CO was then introduced into 
the anode mixture for a further 1h. Beyond t=2h the CO 
concentration was alternately set to zero for 10 min and to suc- 
cessively increased values for 10 min, eventually set to zero until 
t = 5h. For Veen = 0.6 V the calculations were then repeated 
with 40%H?2/0.1%O2/N2 and 40%H2/4%O2/N2. The results are 
demonstrated in Fig. 11, the left-hand figure at three differ- 
ent cell voltages without bleeding and the right-hand figure 
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at two different O2 concentrations (in the anode channel) for 
Vece = 0.6 V. From these plots we can make the following 
observations. 


(1) The cell has a capacity to recover quite well by simply 
reversing to 40% Hz, even at a low cell voltage, as observed 
in the experiments in [29,50,51]. 

On these relatively short time scales O2 bleeding even at 
4% does not significantly counter the poisoning effect when 
the CO concentration is high (particularly above 100 ppm), 
in agreement with the experimental results in [29,50,51]. 
Fig. 12 shows the evolutions of surface coverages eval- 
uated at the membrane/ACL interface for the two anode 
oxygen concentrations; the CO coverages in the two cases 
are increasingly closer as the CO concentration is increased. 
Although not a proof, this latter result provides an indication 
that, as is known from experiment, there is an upper limit to 
the ameliorating effect of bleeding. 
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Fig. 13. The left hand figure shows polarization curves (from potential sweeps of 1-0 V in 80 min) with 25 ppm CO in 40%H2/N2/02, comparing bleeding with 
non-bleeding at two channel temperatures. The remaining parameter values are given in Tables 5—7. The right-hand figure shows the evolutions of surface coverages 


evaluated at the interface between the ACL and membrane, x = x3, in both cases. 
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A more striking feature of the transient operation is seen 
when the channel temperatures are varied. The left-hand side 
of Fig. 13 compares bleeding with non-bleeding at two chan- 
nel temperatures (same in both channels). These polarization 
curves are generated from simulated potential sweeps (1 V to 
0 V in 80min) with 25 ppm CO in 40%H2/N2/O2. The perfor- 
mance gain from bleeding is markedly greater at 70 °C, the lower 
temperature. Notice the shape of the two curves representing 
non-bleeding: the better performance in the upper part of the 
curve for the lower temperature is almost certainly influenced 
by the temperature dependence of the adsorption and desorption 
rate constants (see the preceding discussion on Fig. 10 and the 
right-hand side of Fig. 9). The same trend is observed in the two 
cases representing bleeding, but with the crossover of the two 
curves happening at a much lower cell voltage. The right-hand 
side of Fig. 13 shows the evolution of CO, H2 and O2 coverages 
for the two cases with bleeding. Recall that a lower temperature 
leads to lower adsorption and desorption rate constants for H2 
and CO with the decrease in the CO adsorption constant domi- 
nating. Moreover, the ratio of O2 desorption to adsorption rate 
constant, Bis is decreased by an order of magnitude at the lower 
temperature. We conjecture that this dependence (and that of the 
protonic conductivity) on temperature influence the result; how- 
ever, system (25)—(27) is simply too non-linear to unravel with 
a single example the exact path that leads to the profiles and 
curves in Fig. 13. Adsorption, desorption, electro-oxidation and 
heterogeneous catalysis also depend non-linearly on the surface 
coverages, concentrations and overpotential, which are in turn 
influenced by other quantities. This point is not pursued fur- 
ther (apart from a brief discussion later) but is worthy of further 
exploration, preferably in combination with experiment. 


4.2. Effects of Pt-Ru alloys 


By incorporating the bi-functional mechanism proposed by 
Watanabe and Motoo [39], in this section we investigate the 
role of Ru in designing ‘CO tolerant’ catalyst alloys. Modified 
versions of this mechanism have previously been employed for 
PEMFC, most notably by Camara et al. [35] and by Enbiack and 
Lindbergh [36]. Both cases are restricted to the kinetic problem 
in the anode, the former including H; (at steady-state) and the lat- 
ter neglecting it (with transient effects included). Here we use the 
dualistic approach in [36], in which it is assumed that CO adsorbs 
preferentially onto the Pt sites and water adsorbs preferentially 
onto the Ru sites, and that the adsorbed OH species are mobile 
enough to reach the Ru sites where CO is oxidized. It is assumed 
therefore that the surface electro-oxidation of CO by OH is 
dependent on the CO coverage on Pt and OH coverage on Ru. It 
is known that H2 electro-oxidation on Ru is roughly two orders 
of magnitude smaller than on Pt [32], and is therefore neglected. 
Adsorption of water occurs at a significantly lower potential on 
pure Ru than on Pt [52], and so adsorption on Pt is not included 
as a simplification. Adsorption of CO on Ru is also considered 
small in this approximation, as in [36]. Though somewhat sim- 
plified, we shall demonstrate that this model captures several 
salient features and is qualitatively similar to a more complex 
model that includes CO adsorption and oxidation on Ru. 


The (bi-functional) mechanism to be modelled (supplement- 
ing (23)) is 


HO adsorption on Ru: H20 + M = OHM + H+ +e, 


CO electro-oxidation by OH on Pt : 
CO-M + OH-M = Ht +e” + CO2 + 2M (34) 


In this dualistic approach the already-defined surface coverages 
are now assumed to be valid on the Pt sites only. The Ru sites 
are assumed to be occupied by OH, which has surface coverage 
qm. or else free, with surface coverage ak 


6co + 64 + Om = 1, ooa + OM = 1 


assuming that there is no O2 present. The equations for the 
evolution of the coverages of CO (on Pt) and OH (on Ru) are 
therefore 


d6co 
P— a, = ICO,ads — YCO,ox ~ COOH, ox (35) 
Ru 
PRu = H)0,ads — YCO—OH,ox (36) 


where the extra terms (arising from (34)) are defined as 


‘4 Ru 
4CO-OH,ox = k, 08". co ef a/2RT QBroco/RT ob2r200g/ RT 


—k_ OR" OmCco,,¢ e` Filta/2RT .—(U1—B)r6c0/RT 


x e7 (1-22r200u/RT (37) 


and 


= Ru 
4H)0,ads = koCgoRe ef Ma/2RT o B2r200g/ RT 


—k_26R4, e7 Fita/2RT (1—2)r265),/ RT (38) 


It is assumed that the adsorption of water on Ru occurs from 
the dissolved phase, so that the right-hand side of Eq. (8) now 
contains an extra term in the ACL. Similarly, the ACL terms 
Sg, and Sg, in Eq. (7) are modified to account for proton and 
electron generation in (34). Therefore, in Table 2 we now have 
(for the ACL) 


RHS of Eq. (8) = — Saq — Sal — @RugH,0,ads; 


So. = — So, = 4F (Gu.ox + 24C0,ox + YCO-OH)+4Ru GH 0, ads 
(39) 


Moreover, the active surface area of Pt, a, is now a function of 
the Pt—Ru ratio, rp—Ru 


a = ATrPt—Ru> aRu = 4T —a 


where ar is the combined active surface area of Pt and Ru and 
aru is the surface area of Ru. 

The kinetics for H2O adsorption are based on a Frumkin 
isotherm, [53], with interaction parameter, r2, and symmetry 
coefficient, 62, both of which were fixed at the values given in 
[36]. The values of the forward and backward rate constants k2 
and k_2 were also based on those in [36] (see Table 6). In order 
to achieve the correct qualitative trends the parameters kų and 
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Fig. 14. A demonstration of effects of changes in the ratio of platinum to ruthenium and in the channel temperatures. For a list of the remaining parameters and 


conditions see Tables 6 and 7. 


k_, were adjusted. The selected value of kų was much larger 
than that used in [36], reflecting the competition between CO 
oxidation and the other mechanisms in (23) and (34), as well as 
the influence of the transport phenomena. Neither of these were 
not considered in [36]. 

Fig. 14 shows the effect of changes in both the ratio of Pt 
to Ru and in the channel temperatures (other parameters as in 
Tables 5-7). These are again simulations of potential sweeps 
from 1 to OV in 80min. The best performance is seen for a 
Pt—Ru ratio of 1—1 and the worst for pure Pt, with the perfor- 
mance for a Pt—Ru ratio of 1—9 intermediate between these two 
case, at both temperatures. The use of Pt—Ru involves a trade- 
off between increased CO electro-oxidation (by the OH species) 
occurring at low potentials on the Ru, and decreased H2 oxida- 
tion on the Ru portion of the catalyst surface, implying that 
there is an optimal ratio of Pt to Ru surface area (or volume) 
at which the best performance on reformate is achieved. This 
result is in agreement with the experiments in [54] (and refer- 
ences therein), where the authors concluded an optimal Pt-Ru 
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ratio of between 1-1 and 1-3. Moreover, it supports the appli- 
cability of the bifunctional mechanism (34). The results shown 
in Fig. 14 suggest that at elevated temperatures the effect of Ru 
on the performance is less visible, matching the conclusion in 
[54]. The reasons for this will include the temperature effect on 
membrane conductivity and on CO adsorption and desorption 
already encountered in Fig. 9. 

Adsorption of CO on Ru was neglected in the model above. 
Assuming that adsorption and oxidation can be modelled with 
Eq. (26), i.e. the following equation for CO on Ru is added to 
(35): 


dogo 2 =k e7 rol RT (Cog, ore 


PRu ~ dt 


Ru 
— beO POR e rc o/ R ry 


Fa 
19 (40) 


—2k co sinh ( 


the results are qualitatively the same if the rate constants 
are adjusted. In Fig. 15 the results at 60°C in Fig. 14 are 
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Fig. 15. A comparison of the models with CO adsorption and oxidation on Ru included (left-hand side) and neglected (right-hand side). 
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compared with a model that includes Eq. (40). The rate con- 
stants were fitted ask, = 2.5 x 1076 mol m~? s~!and ky = 1x 
1075 mol m~? s~!(the other constants as in Table 6), which are 
much larger than the values in Table 6 since the H20 now com- 
petes with CO for adsorption onto Ru. 


5. Conclusions 


CO poisoning and Oz bleeding in PEMFC are largely influ- 
enced by the (rates of) the electrochemical reactions in the anode. 
However, neglecting other important effects reduces the number 
of realistic settings that we can hope to simulate with modelling; 
these effects can influence and also be influenced by the kinet- 
ics. This point is amply illustrated in the examples provided, 
such as Figs. 7 and 8 and the associated discussion; exclud- 
ing liquid-water transport and membrane hydration can result in 
conclusions that hold true only in the most idealized of settings. 
In the general case PEMFC operation is the conflation of several 
competing processes, and if a model is to be used as a predictive 
or diagnostic tool inclusion of these fundamental phenomena is 
likely to be necessary. 

Equally as important for automobile applications (in partic- 
ular) is that steady-state modelling of CO poisoning provides 
incomplete, and perhaps misleading information. In, for exam- 
ple, the right-hand side of Fig. 9 and the left-hand side of Fig. 13 
we see that evaluating the response of the system on a long time 
scale can mask shorter term features that are of potential signifi- 
cance for real events that can last for seconds, minutes or even a 
few hours. The results contained in these examples, demonstrat- 
ing the effects of temperature on the transient response, could 
have far-reaching implications if verified by experiment. 

The parameter space of possible investigation is extremely 
high dimensional, reflecting the many physical properties and 
external influences that can be important in the phenomenology. 
In the current model Tables 5-7 list about 75, many of which take 
a value specific to a particular physical configuration or mate- 
rial. Enumeration of all possibilities is therefore impracticable. 
The primary objective has been to demonstrate the flexibility in 
the framework and its ability to capture, predict and possibly 
explain the complex interactions observed in experiment. Nat- 
ural techniques to investigate the short- and long-term effects 
of CO poisoning and O2 bleeding include cyclic voltamme- 
try, as in [49], and the pulsing experiments found in [50,51]. 
In principle, such experiments are less time-consuming than 
steady-state alternatives, and are more applicable to PEM fuel 
cell operation in many if not most applications. The modelling 
tool we have outlined in this paper is an ideal accompaniment 
for such experiments, to help reduce the number (and therefore 
cost) of laboratory investigations, and to help understand the 
results. 

Of course, there are several stages in reaching such an out- 
come. There remains some uncertainty in the (electro) chemical 
rate constants (29). A great deal more transient data is required 
to estimate or fit these parameters in as wide a range of param- 
eter space as possible, including investigations into the effects 
of sweep rate, temperature, anode pressure and water activity. 
Characterization of the electrochemistry for Pt alloys (other than 


Pt-Ru) currently in use would also be of extreme benefit. These 
avenues, along with a confirmation of the widely occurring 
behaviour demonstrated in Figs. 9 and 13, represent practical 
complements to the present work. 
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